function FVlinearadvectionFOUpseudo1D 
% File name: FVlinearadvectionFOUpseudo1D.m 
% Description:  Solves the PDE, dU/dt + div(H) = 0 
% using the FOU finite volume scheme on a uniform Cartesian mesh inclined 
% at theta degrees to the x-axis.  This is a pseudo-1D problem.  There is 
% a single row of NI rectangular cells in the i direction with side lengths 
% di and dj.  Flow velocity v has been fixed to ensure that there is no flow 
% in the j direction. 
% 
% This equation corresponds to the finite volume form: 
% doubleintegral(dU/dt dA) + lineintegral(H.n ds) = 0. 
% 
% In this case the equation is the linear advection equation so:  
% velocity: v = vx i + vy j where vx, vy are constant. 
% flux density: H = vU = vx U i + vy U j. 
% 
% Cell areas and side vectors are the same for each cell. 
% 
% The program is written in a structured way using subfunctions so that it 
% can be modified easily to pseudo-1D or 2D casees on non-Cartesian 
% structured meshes. 
% 
% Initial conditions:  Gaussian along 1D line of cell centres. 
% 
% Boundary conditions:  Left: zero gradient. 
% 
% subfunctions: freadmesh, fcellarea, fsidevectors, finitialu, 
%               finterfaceflux, fghostcells, fcalcdt, fplotresults,  
%               fdrawmesh 
% 
% Verification: Looks correct - numerical and exact solutions agree 
% and the profile has travelled the correct distance.  Proper verification 
% requires a pen and paper calculation on a small mesh. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% 
clc; clf; clear all; 
runtime=10.0; % runtime in seconds 
t=0;          % current time 
timelevel=0;  % current time level 
[x,y,xcen,ycen,solid,v]=freadmesh; % gets coordinates of the 1D computational mesh  
              % and solid flags (excluding ghost cells) and flow velocity 
              % Note that there are no solid interior cells in 1D. 
disp('mesh read in') 
[m,n]=size(x);  
NI=m-1; % number of computational cells in i direction. 
A=fcellarea(x,y);   % computes constant cell area 
disp('calculated cell area') 
S=fsidevectors(x,y); % compute and store cell side vectors                                          
disp('calculated cell side vectors') 
u0=finitialu(xcen,ycen,t,v,S); % Put initial cell centre values of u in a NIx1 array 
uinitial=u0; % store initial profile for plotting 
disp('inserted initial conditions') 
% Append extra cells around arrays to store any ghost values. 
u0=finsertghostcells(u0);  % u0 is now (NI+2)x1 
u1=zeros(size(u0));  % (NI+2)x1 array for u values at next time level 
%% 
disp('time marching starts')    
while(t<runtime)    
 timelevel=timelevel+1   
 u0=fbcs(u0);  % Implement boundary conditions using ghost cells.  
               % u0 is (NI+2)x1 and each computational cell is  
               % indexed i=2 to NI+1. 
 dt=fcalcdt(A,S,v); % Finds stable time step for each iteration. 
 t=t+dt  % update time   
   for i=2:NI+1  
       IH=finterfaceflux(v,u0,i); % gets the left and right interface fluxes for cell i 
       IHright=[IH(1,1),IH(1,2)]; % interface flux vector right side  
       IHleft=[IH(2,1),IH(2,2)];  % interface flux vector for left side 
       % 
       sright=[S(1,1),S(1,2)]; % side vector for right side 
       sleft=[S(2,1),S(2,2)];  % side vector for left side 
       % FV scheme 
       % compute total flux out of cell i 
       totalfluxout=dot(IHright,sright)+dot(IHleft,sleft); 
        
       % totalfluxout=vx*u0(i,1)*dy+vx*u0(i-1,1)*(-dy); WORKS!!!!!!! 
       u1(i,1)=u0(i,1)-(dt/A)*totalfluxout;  
   end % of i loop 
 %        
   u0=u1;  % update u values 
end  % of while loop 
uexact=finitialu(xcen,ycen,t,v,S); % find exact solution 
t  % print time 
subplot(2,1,1),fdrawmesh(x,y,solid)   
%fdisplay(u0); % print results for a small mesh as a matrix 
subplot(2,1,2), fplotresults(xcen,ycen,u0,uinitial,uexact) % plot results 
disp('program ended at time t see plots') 
end % FVlinearadvectionFOUpseudo1D 
%%%%%%%%%%%%%%%%%%%%%%%%%%% subfunctions  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [x,y,xcen,ycen,solid,v]=freadmesh 
% Verification:  looks correct using fdrawmesh 
% The mesh is structured and has NI cells in the i direction and  
% 1 cell in the j direction. 
% The x and y coordinates of the lower left hand corner of cell (i,j) are 
% held in arrays x and y respectively which are both (NI+1)x2.  In 
% this way the 4 vertices of cell (i,j) are (x(i),y(j)), (x(i+1),y(j)), 
% (x(i+1),y(j+1)) and (x(i),y(j+1)).  solid is an NI by 1 array which  
% flags solid cells.  If cell (i,j) is solid then solid(i,j)=1 otherwise 
% solid(i,j)=0. 
% x(1,1)=0=y(1,1). 
% The mesh is created as a single row of uniform cells in the i direction 
% which is a straight line inclined at theta degrees to the x axis.   
% This problem is pseudo-1D. 
% v=[vx,vy]: vy component is chosen so that there is no flow in  
% the j direction. 
%  
NI=100;  % number of cells in the i direction 
x=zeros(NI+1,2); % allocate correct sized array for cell vertices 
y=zeros(NI+1,2); % allocate correct sized array for cell vertices 
xcen=zeros(NI,1); % allocate correct sized array for cell centres 
ycen=zeros(NI,1); % allocate correct sized array for cell centres 
solid=zeros(NI,1); % allocate correct sized array for solid cell flags 
di=1.0; % cell length in i direction 
dj=2.0; % cell length in j direction (arbitrary) 
theta=pi/4;  % inclination angle (anticlockwise) 
m=tan(theta); % slope of line 
vx=sqrt(2);   % velocity component in x direction 
vy=m*vx; % correct velocity component in y direction 
% note that for theta=pi/4 the above ensure that flow is only in the 
% i direction and the flow speed is 2 m/s. 
v=[vx,vy]; % flow velocity vector. 
dx=di*cos(theta);  % mesh increment in x direction 
dy=di*sin(theta);  % mesh increment in y direction 
x(1,1)=0; % first row starting point for x 
y(1,1)=0; % first row starting point for x 
x(1,2)=-dj*sin(theta); % second row starting point for x 
y(1,2)=dj*cos(theta);  % second row starting point for y 
% create uniform mesh 
for i=2:NI+1 
      x(i,1)=x(1,1)+(i-1)*dx; 
      y(i,1)=y(1,1)+(i-1)*dy; 
      x(i,2)=x(1,2)+(i-1)*dx; 
      y(i,2)=y(1,2)+(i-1)*dy; 
end % of i loop 
%  x 
%  y 
% find cell centres by averaging coordinates of vertices  
% (this works for a general structured mesh) 
for i=1:NI 
     xcen(i,1)=(x(i,1)+x(i+1,1)+x(i+1,2)+x(i,2))/4; 
     ycen(i,1)=(y(i,1)+y(i+1,1)+y(i+1,2)+y(i,2))/4; 
end % of i loop 
end % freadmesh 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
function A=fcellarea(x,y)   
% Verification: works but needs to be generalised for non-rectangular cells 
% In a uniform Cartesian mesh cell area = di dj  
di=sqrt((x(2,1)-x(1,1))^2+(y(2,1)-y(1,1))^2); 
dj=sqrt((x(1,2)-x(1,1))^2+(y(1,2)-y(1,1))^2);    
A=di*dj;   
end % fcellarea 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
function S=fsidevectors(x,y) 
% Verification: not verified 
% For each cell: calculates right and left side vectors. 
% 
S=zeros(2,2); 
% 
rightside=[x(2,2)-x(2,1), y(2,2)-y(2,1)]; 
% Constant right side vectors 
S(1,1)=rightside(2); 
S(1,2)=-rightside(1); 
% Constant left side vectors 
leftside=[x(1,1)-x(1,2), y(1,1)-y(1,2)]; 
% Constant right side and unit normal vectors 
S(2,1)=leftside(2); 
S(2,2)=-leftside(1); 
end % fsidevectors 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
function u=finitialu(xcen,ycen,t,v,S)  
% Verification:  verified 
% Inserts u values at time t.  
[NI,NJ]=size(xcen); 
u=zeros(NI,1); 
rightside=[S(1,1),S(1,2)]; 
length=sqrt(dot(rightside,rightside)); % length of side vector 
rightunitnormal=rightside/length; 
speed=dot(v,rightunitnormal);  % flow speed in i direction 
% Initial 1D Gaussian function of d where d is the distance along the line  
% connecting cell centres from (xcen(1,1),ycen(1,1)). Gaussian is based at d/2; 
dmax=sqrt(((xcen(NI,1)-xcen(1,1))^2+(ycen(NI,1)-ycen(1,1))^2)); % max d value of domain 
dc=dmax/3;  % approx mesh centre distance 
% 
cutoff=dmax/4 % cut off radius for Gaussian 
for i=1:NI 
        d=sqrt((xcen(i,1)-xcen(1,1))^2+(ycen(i,1)-ycen(1,1))^2)-speed*t; 
        if (abs(d-dc)<cutoff) 
           u(i,1)=exp(-0.01*(d-dc)^2);  % Gaussian 
        else 
           u(i,1)=0; 
        end 
end % of i loop 
end % finitialu 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function IH=finterfaceflux(v,u0,i) 
% Verification:  not verified 
% Calculates right and left interface fluxes for each cell. 
% Flux H depends on U, i.e. H = H(U).  
% For the 1D linear advection equation, H(U) = vx U i + 0 j. 
%  
% There is considerable choice for interface flux estimation: different  
% choices give different FV schemes. 
%  
%% The following FOU scheme simply takes: 
%  IHright=H(u0(i,1)) = v u0(i,1) 
%  IHleft=H(u0(i-1,1)) = v u0(i-1,1) 
%  This works (is an upwind scheme) for vx > 0. 
% 
IH=zeros(2,2); % array to store left and right interface flux vectors 
% 
% Right 
IH(1,1)=v(1)*u0(i,1);  % i component of right interface flux 
IH(1,2)=v(2)*u0(i,1);  % j component of right interface flux 
% Left 
IH(2,1)=v(1)*u0(i-1,1);  % i component of left interface flux 
IH(2,2)=v(2)*u0(i-1,1);  % j component of left interface flux 
end % fIflux 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
function outarray=finsertghostcells(inarray) 
% Verification:  not verified 
% Assumes that ghost cells are needed at either end of the domain. 
% array is embedded into (m+2)x1 array of zeros.  Hence computational  
% indices go from i=2 to m+1.  In the FOU scheme a ghost cell 
% is needed only at the left end. 
[m,n]=size(inarray); 
dummy=zeros(m+2,1); 
dummy(2:m+1,1)=inarray; 
outarray=dummy; 
end % finsertghostcells 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
% function outarray=fremoveghostcells(inarray) 
% Verification:  not verified 
% % Removes ghost cells so that a mx1 array becomes (m-2)x1 
% [m,n]=size(inarray); 
% outarray=inarray(2:m-1,1); 
% end % fremoveghostcells 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function u=fbcs(u) 
% Verification:  verified 
% Implements boundary conditions using ghost cells index by i=1, i=NI+2. 
[m,n]=size(u); 
NI=m-2; 
% Left: zero gradient 
u(1,1)=u(2,1); 
% Right: not required for FOU scheme 
end % bcs 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function dt=fcalcdt(A,S,v) 
% Verification:  not verified 
% Finds allowable time step dt using heuristic formula. 
% This formula must be generalised for a non Cartesian mesh. 
rightside=[S(1,1),S(1,2)]; % right side vector 
F=0.95; % tuning factor 
dt=A/abs(dot(v,rightside)); % heuristic formula 
dt=F*dt; 
end  % fcalcdt 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function fdrawmesh(x,y,solid) 
% Verification:  not verified 
% Description:  Draws a structured 2D finite volume mesh and fills in any 
% solid cells.  
% Date structures: 
% The mesh has NI cells in the i direction and NJ cells in the j direction. 
% The x and y coordinates of the lower left hand corner of cell (i,j) are 
% held in arrays x and y respectively which are both NI+1 by NJ+1.  In 
% this way the 4 vertices of cell (i,j) are (x(i),y(j)), (x(i+1),y(j)), 
% (x(i+1),y(j+1)) and (x(i),y(j+1)).  solid is an NI by NJ array which  
% flags solid cells.  If cell (i,j) is solid then solid(i,j)=1 otherwise 
% solid(i,j)=0. 
%  
[m,n]=size(x); 
NI=m-1; % number of cells in i direction 
NJ=1; % number of cells in j direction 
% 
% Plot the mesh 
hold on   % keeps all plots on the same axes 
% draw lines in the i direction 
for i=1:NI+1 
    plot(x(i,:),y(i,:),'k') 
end 
% draw lines in the j direction 
for j=1:NJ+1 
    plot(x(:,j),y(:,j),'k')    
end 
title('computational mesh') 
xlabel('x') 
ylabel('y') 
% Fill in solid cells 
for i=1:NI 
    for j=1:NJ 
        if (solid(i,j)==1) 
            solidx=[x(i,j),x(i+1,j),x(i+1,j+1),x(i,j+1),x(i,j)]; 
            solidy=[y(i,j),y(i+1,j),y(i+1,j+1),y(i,j+1),y(i,j)]; 
            fill(solidx,solidy,'k')       
        end 
    end % of j loop 
end % of i loop 
end % fdrawmesh 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function fplotresults(xcen,ycen,u0,uinitial,uexact) 
% Verification:  verified 
% Dispays results along the mesh (i direction) 
[NI,NJ]=size(xcen); 
u0=u0(2:NI+1,1); % extract computational values 
% find distances along the cell centre line 
for i=1:NI 
    d(i)=sqrt((xcen(i,1)-xcen(1,1))^2+(ycen(i,1)-ycen(1,1))^2); 
end    
plot(d,uinitial,'k--',d,uexact,d,u0,'k.') 
title('Graphs of initial profile (--) and numerical(.) and exact solutions') 
xlabel('distance [m] along cell centre line of the pseudo-1D mesh') 
ylabel('U [kg/m]') 
end % fplotresults 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 

